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ABSTRACT 

Three-dimensional (3 D) hydrodynamic simulations of shell oxygen burning 
(jMeakin ArnettI l2007bl ) exhibit bursty, recurrent fluctuations in turbulent kinetic 
energy. These are shown to be due to a general instability of the convective cell, re- 
quiring only a localized source of heating or cooling. Such fluctuations are shown to be 
suppressed in simulations of stellar evolution which use mixing-length theory (MLT). 

Quantitat ivelv similar behavior occurs in the model of a convective roll (cell) of 



Lorena (jl963l ). which is known to have a strange attractor that gives rise to chaotic 
fluctuations in time of velocity and, as we show, luminosity. Study of simulations 
suggests that the behavior of a Lorenz convective roll may resemble that of a cell in 
convective flow. We examine some implications of this simplest approximation, and 
suggest paths for improvement. 

Using the Lorenz model as representative of a convective cell, a multiple-cell model 
of a convective layer gives total lumino sitv fluctuations wh ich are suggestive of irregular 
variables (red giants and supergiants (ISchwarzschildlll975l)). and of the long secondary- 
perio d feature in semiregular AGB variables ( Stother J 1201(1 : IWood. Olivier fc Kawalerl 
20041 ). This "r-mechanism" is a new source for stellar variability, which is inher- 
ently non- linear (unseen in linear stability analysis), and one closely related to in- 
termittency in turbulence. I t was already implicit in the 3D global simulations of 
Woodward Porter, fc Jacobsl (|2003l ). T his fluctuating behavior is seen in extended 2D 
simulations of CNeOSi burning shells (jArnett &: MeakinI l2011bl ). and may cause in- 
stability which leads to eruptions in progenitors of core collapse supernovae prior to 
collapse. 
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1. Introduction 

Three-dimensional fluid dynamic simula- 
tions of turbulent convection in an oxygen- 
burning shell of a presupernova star show 
bursty fluctuations which are not seen in one- 
dimensional stellar evolutionary calculations 
(which use va rious versions o f mixin g-length 
theory, MLT, lBohm-Viteris3 (jigssl l). This 
paper explores the underlying physics of this 
new phenomena. 



Since the formulation of MLT (jBohm-Vitense 

19581 ). there have been a number of significant 
developments in the theoretical understand- 
ing of turbulent convective flow. 

Fir st, Kolmogorov ( 19621 ) and Qbukhov 
(|l962l ) developed the modern version of the 
turbulent cascade, and published in journals 
easil y accessible in the west; the original the- 
ory ( Kolmogorov 194ll ) was not used in MLT 
although it pre-dated it. This explicit expres- 
sion for dissipation of turbulent velocitie^. 



^turb 



(1) 



where Urms is the root-mean-square of the 
turbulent velocity and is the dissipation 
length. It is found both experimentally and 
numerically that w ^cz-, where tcz is 
the depth of the convective zone. Simula- 
tions for low-Mach number flow show that 
the average of this dissipation over the con- 
vective zone closely compensates for the cor- 
responding average of the buo yant power 
(jArnett. Meakin. h You^^ l2009l ). This ad- 
ditional constraint allows an alternative to 
present practice: fixing the free parameter 
(e.g., the mixing length factor a) directly by 
terrestrial experiments and numerical simu- 
lations whic h deal with the process of turbu - 
lence itself ( Arnett. Meakin. fc Young 20ld ). 
instead of calibrating it from complex astro- 



^Thisis essen tially the "four-fifths law of Kolmogorov" , 
iFrischI (|l995l ). p. 76. 



nomical systems (stellar atmospheres) as is 
now done. 

Second, there has been a considerable de- 
velopment in understanding the nature of 
chaotic be h avior in nonlinear systems; see 
Cvitanovid (|l989l ) for a review and r eprints 



of ori gi nal papers, and Fris ch^ ('1995II: iGleickl 



1987 ) ; [Thompson &: Steward. (19861 ) . iLorenz 
19631 ) presented a simplified solution to 



the Rayleigh pr o blem of thermal convection 
( Chandrasekhai] Il96ll ) which captured the 
seed of chaos in the Lorenz attractor, and 
contains a representation of the fluctuating 
aspect of turbulence not present in MLT. 
This advance was allowed by the steady in- 
crease in computer power and accessibility, 
which lead to the exploration of solutions for 
simple sys tems of nori li near d ifferential equa- 
tions (see Cvitanovic ( 19891 ) and references 
therein) 



It became clear that the Landau 



picture (jLandaul 1 1944 ) of the approach to 
turbulence was incorre ct both theoretically 



( Ruelle &: TakensI Il97lh a nd experimentally 
(|Libchaber Mauerl 1198^ 1. A striking fea- 



ture of these advances has been the use of 
simple mathematical models, which capture 
the essence of chaos in a model with much re- 
duced dimensionality compared to the physi- 
cal system of interest. 

Third, it has become possible to simulate 
turbulence on computers. This realizes the 



vision of John von Neumann (jvon Neumann 
19481 ). in which numerical solutions of the 



Navier-Stokes equations by computer are used 
to inform mathematical analysis of turbu- 
lence. In this paper we will follow this idea 
of von Neumann, in the style which proved 
successful for chaos studies: building simple 
mathematical models of a more complex phys- 
ical system (in this case, the numerical simula- 
tions of turbulent convection) . This approach 
should lead to algorithms suitable for imple- 
mentation into stellar evolution codes, which, 
unlike MLT, are (1) based upon solutions to 
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fluid dynamics equations, (2) non-local, (3) 
time-dependent, and (4) falsifiable by terres- 
trial experiment and future improved simula- 
tions. 

Our particular example is a set of simula- 
tions of o xygen burning in a shell of a star 
of 23Mq (jMeakin fc Arnetd l2007bl l. This is 
of astronomical interest in its own right as 
a model for a supernova progenitor, but also 
happens to represent a relatively simple and 
computationally efficient case, and has gen- 
eral implications for the convection process in 
all stars. 

Three-dimensional hydro dynamic simula- 
tions of shell oxygen burning exhibit bursty, 
recur rent fluctuations in turbule nt kinetic en- 
ergy (jMeakin fc ArnettI (|2007bl ) and below). 
The reason for this behavior has not been ex- 
plained theoretically. These simulations show 
a damping, and eventual cessation, of tur- 
bulent motion if we artificially turn off the 
nucle ar burning ([Arnett. Meakin. &: Yound 
2009 ). Further investigation (jMeakin fc Arnett 
2003 ) shows that nearly identical pulsations 
are obtained with a volumetric energy gener- 
ation rate which is constant in time, so that 
the cause of the pulsation is independent of 
any temperature or composition dependence 
in the oxygen burning rate. Localized heating 
is necessary to drive the convection; even with 
this time-independent rate of heating, pulses 
in the turbulent kinetic energy still occur. 

Such behavior is fundamentally differ- 
ent from traditional nuclear-energized pulsa- 



tions dealt with in the literature 
£-mec hanism. 



Ledouxl (|l94ll ): lUnno. et al 



^e-g-^ the 



and is a consequence of time-dependent 
turbulent convection (it might be called a "r- 
mechanism" , with r standing for turbulence) . 
It appea rs to be relevant to all stellar coii - 
vection. Woodward. Porter. Sz Jacobs! (j2003l ) 
found, in a very different context, that non- 
linear interaction of the largest modes excited 



pulsations of a red-giant envelop^, which is 
another example of the r-mechanism. 

In Section 2 we examine the the physics 
context of the turbulence, including impli- 
cations of subgrid and turbulent dissipation 
for the implicit large eddy simulations (ILES) 
upon which our analysis is based, and the ef- 
fect of the convective Mach number on the 
nature of the flow. In Section 3 we review 
the 3D numerical results of shell oxygen burn- 
ing which are relevant to the theory. In Sec- 
tion 4 we prese nt the results of the classical 
Lorenz model (Lorenz 19631 ) for conditions 
similar to those in Section 3. In Section 5 
we consider implications of turbulent inter- 
mittency on stellar variability, and provide a 
model light curve from this effect alone. Sec- 
tion 6 summarizes the results. The appendix 
gives a short derivation of the Lorenz model. 

2. Physics Context. 

In this section we summarize concepts 
which are needed for the interpretation of 
later results. 

2.1. Subgrid Damping and Kolmogorov 

Approximation of partial differential equa- 
tions by discrete methods inevitably leads to 
a loss of information at scales smaller than the 
grid size. A single element in space is approxi- 
mated as a homogeneous entit}0; this is equiv- 
alent to complete mixing at this scale, at each 
time step, of mass, momentum, and energy. 
The loss of information that occurs with this 
mixing co rresponds to an increase in entropy 
(j Shannon! 1 1948) . the mixing of momentum is 



K-mechanism, which depends upon variations in 
opacity, is not required to drive such pulsations. 
^Actually most modern simulations (ours included) use 
higher order methods which make some further as- 
sumptions regarding the behavior of variables inside a 
zone. This complicates but does not change the ar- 
gument; it is still true that information is lost at the 
zone level. 
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equivalent to the action of viscosity, and the 
mixing of internal energy corresponds to th e 
transport of heat ([Landau &: Lifshitz (Il959l l. 
§15 and §49). 

In 3D flow, turbulent energy will cascade 
from large scal es to small, a, t a rat e set by the 
largest scales (iKolmog orov 1941 ). At suffi- 
ciently small scales, microscopic processes ho- 
mogenize the flow and dissipate the kinetic 
energy. Thus, there is a deep connection 
between the turbulent cascade and sub-grid 
scale mixing. 

(|200d ) 



Svtine. et al. 



have demonstrated 
that the piece- wise parabolic method (PPM, 
which we use), based on the Euler equation 
(which has no explicit viscosity), converges to 
the same limit as methods based on Navier- 
Stokes equation (which do have explicit vis- 
cosity), as the grid is refined to smaller zones 
and smaller effective viscosity (the relevant 
limit for astrophysics). The subgrid scale dis- 
sipation for monoto nicity preserving hydro- 
dyna mic algorithms ( Boris! 2007 : Woodward 
20071 ). which is implicit in these methods, au- 
tomatically gives a reasonable treatment of 
the turbulent cascade down to the grid scale. 
We use this implicit sub-grid dissipation in 
our large eddy simulation (ILES); this is the 
most computationally efficient way to deal 
with turbulent systems with a large range of 
scales. The largest scales, which set the rate 
of cascade and contain most of the energy are 
explicitly calculated, while the sub- grid scales 
are dissipated in a way consistent with the 
Kolmogorov cascade. 



Woodward et al.l (|2006l ) have presented a 



refinement of the PPM algorithm which has 
improved behavior at the smaller resolvable 
scales. We have not yet implemented this 
modification. The theoretical approach used 
here involves integrated properties of convec- 
tive cells; we find by direct resolution stud- 
ies that these properties are well estimated 
even with surprisingly modest resolution, be- 



cause they are determined primarily by the 
largest scales in the convective region. It ap- 
pears that the ILES simulations are adequate 
for the present analysis. 

Arnett. Meakin. &: Young (I2OO9I ) have shown 



that the numerical damping at sub-grid scales 
in our ILES simulations is quantitatively con- 
sistent with the introduction analytically of 
the Kolmogorov cascade into the theoretical 
discussion. The turbulent velocity field was 
found to be dominated by two components: 

1. a non- isotropic flow of the largest scale 
modes in the convection zone, which 
is coupled to the fluctuations. This 
has aspects of a "coherent struc ture" 
(jHolmes. Lumlev. fc Berkoodll996l ). The 
largest scales are unstable toward break- 
up, but are least affected by dissipation, 
and in this sense the most laminar. 

2. a more isotropic, homogeneous turbu- 
lent flow which carries the kinetic en- 
ergy via the turbulent cascade to scales 
small enough for d issipation to occur 
( Kolmogorov Because of the 
vast size difference, the small scales are 
weakly coupled to the largest scales, 
which determine the rate at which en- 
ergy flows through the cascade. 

If we approximate the non-isotropic compo- 
nent of the flow (the largest scale of con- 
vection) with that described by the Lorenz 
model, this interpretation captures the oxygen- 
burning fluctuations. 

2.2. Types of Flow 

There are two limiting cases for convective 
flow, depending upon the convective Mach 
number Mconv (the ratio of the fluid speed 
to the local sound speed); these are usually 
termed the "incompressible" {M.conv < 1) 
and "compressible" {Aiconv ~ 1) regimes. 
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Stars are stratified in density, so that the no- 
tion of "incompressibhty" is misleading. We 
will use the term "low-Mach number flow" in 
place of "incompressible" when we mean flows 
in which acoustic radiation is small, but may 
be compressed due to stratification. 

For turbulent motion, the pressure pertur- 
bation P' is related to the convective Mach 
number by P'/P ~ pulms/P ~ -^coni,- 
Sound waves outstrip fluid motion, so that 
pressure differences quickly become small, ex- 
cept possibly for a static background strat- 
ification. Most of the historical research 
on convection (e .g., t he Benard problem, 
Chandrasekhail ( 1961 ): Landau &: Lifshitz 
( I959I )) is done in this limit. 

Using the Reynolds decomposition, ip = 
+ if' , with horizontal averaging {(p) = 999 
and ((/?') = 0, mass conservation for a steady 
state flow can be written 



(pu) = (pouo) + {p'u') = 0. 



(2) 



The Navier-Stokes equation, using mass 
conservation, is 

dtpvi + V-puu = -VP + pg + z^V^u, (3) 

where puu is the Reynolds stress tensor. Mass 
conservation implies a convenient identity in- 
volving total CO- moving derivatives. 



DtpvL = pDtU. 



(4) 



Taking the dot product of the velocity vec- 
tor u with the Navier-Stokes equation, gives 
a kinetic energy equation, 

Dtpu ■ u/2 = -u • VP+pu ■ g+z/u-V^u. (5) 

If on average the system is in a steady state, 
the time derivative must integrate to zero over 
the convective region, and the mass conserva- 
tion law implies that the total buoyancy power 
term is zero, {pu • g) = (assuming constant 
g on horizontal averaging) , and therefore does 



not contribute to the production of kinetic en- 
ergy anywhere in the flow. The only other 
term, which remains to balance the viscous 
dissipation of the kinetic energy, is the pres- 
sure term 



u • VP = uo • VPo + (u' • VP') 
which may be rewritten as 



(6) 



u • VP = uo • VPo + (V • (P'u')) - (P'V • u'). 

(7) 

The divergence term vanishes upon integra- 
tion over the volume. Using hydrostatic equi- 
librium for the background state, VPq = pog? 
and mass conservation, PqUq = —{p'u'), 

(-U • VP) = (p'u' • g) + (P'V • u') (8) 

When the Mach number is small, the second 
term the right hand expression is nearly zero 
because V • u' ^ and the turbulent pres- 
sure is negligible. In this limit the kinetic en- 
ergy production is best understood as due to 
the remaining buoyancy power term (p'u' • g) , 
which is directly related to the enthalpy flux 
(|Arnett. Meakin. k YoTm3l2009l ^. 

When the Mach number is no longer small, 
the second term on the RHS of Eq. [8] in- 
creases in importance: both the divergence of 
the fluctuating velocity field and the pressure 
perturbation begin to play a role. The ve- 
locity field changes character; it is no longer 
dominated by rotational flow, but develops 
an irrotational component (V • u' 7^ 0). The 
flow becomes diverging (consider the extreme 
limit of a point explosion which is pure di- 
vergence). Also, the ram "pressure" (a tensor 
puu) is not negligible and must be included in 
the momentum equation ("hydrostatic equi- 
librium"). Sound wave g eneration increases 
rapidl y as Adconv 1 (jLandau &: Lifshitz 
( 19591 ). §75). The compressible limit is 
M conv — 1- Shock formation is the most 
startling change in the flow character. 
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Which of the M conv hmits is relevant for 
astrophysics? Both are. Almost all the mat- 
ter in stellar convection zones, during almost 
all evolution, is in the limit of low-Mach 
number flow, as are our turbulence simula- 
tions. The exceptions are important: (1) 
explosions, such as supernovae and novae, 
(2) vigorous thermonuclear flashes, (3) vig- 
orous pulsations, especially radial ones, and 
(4) the sub-photospheric layers of stellar sur- 
face convection zones, which are strongly non- 
adiabatic, to name a few. 

3. The Oxygen Shell Simulation 

Figure [J illustrates the behavior of two im- 
portant integral quantities, the total turbu- 
lent kinetic energy and the total buoyancy 
power, in the oxygen-bur ning shell simula- 
tions for a 23 Mq star (jMeakin &: Arnett 
2007al l. The flow has a low-Mach number 
{Mconv ^ 0.001), although the numerical 
simulations use the equations for fully com- 
pressible fluid flow, and would have correctly 
treated high-Mach number flows. Also shown 
are the same quantities for the Lorenz model, 
which is discussed in Section SI 

At any instant in time, the total convec- 
tive kinetic energy is ^Mczv^ms^ where Mcz 
is the mass of the convective zone and Vrms is 
the rms velocit}0, while the kinetic energy in 
the isotropic part of the turbulent flow field 
is Eturb = ^MczUt, where ut is the turbulent 
velocity, which we define as the isotropic part 
of the turbulent flowj- The reader is warned 



''More precisely, v is the solenoidal velocity in the con- 
vective region, with overall translational velocity re- 
moved, so that Vrms IS its root-mean-square value. 
^This follows I Arnett. Meakin. fc Yound (|2009 l'). §2.4, in 
which the transverse velocities are used to estimate the 
isotropic component of the vertical velocity. Thus, we 

where 



use Ut — - 



4 ^vms 



the velocity components perpendicular to the direc- 
tion of the gravity vector, i.e., the tangential veloci- 
ties. The vertical velocity also contains a significant 
contribution (liirma) from the non- isotropic flow of the 
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Fig. 1. — Comparison of Convective Kinetic 
Energy Fluctuations in the 3D simulations 
and in the Lorenz model, starting from zero 
velocity. The Lorenz model is labeled with 
its Prandtl number a and its Rayleigh num- 
ber r. Turbulent velocity squared (solid) and 
buoyancy power to the 2/3 power (dotted) 
are plotted versus time. The axes where cho- 
sen so that the same number of peaks would 
be shown, and the average kinetic energy be 
the same. The curves are surprisingly simi- 
lar, even though (1) the Lorenz model used 
only a single mode, (2) a and r were not ad- 
justed to give a fit, (3) the Lorenz model im- 
poses thermal balance while thermal imbal- 
ance adds to the background slope in the sim- 
ulation (see Fig. [2]), and (4) the simulations 
have ~ 10^ more degrees of freedom than the 
Lorenz model. 
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that division of the flow into "turbulent" and 
"large scale" flow is useful but an oversimpli- 
fication, so that the exact relative values of 
these two kinetic energies depend upon the 



order unity ( 


Arnett. Meakin. & Youne 


200d; 


Meakin & Arnett 


2011). Consequentlv the 



precise distinction does not change the quali- 
tative picture. 

The buoyancy power is the rate at which 
kinetic energy per unit mass is increased by 
buoyant acceleration. If it is integrated over 
the space containing the convection zone, we 
have Qint = Jczl^^^ where q = —gu'p'/po is 
the buoyant acceleration times the turbulent 
velocity; has units of velocity cubed. 

In Figure [1] the simulations show a phase 
lag of about 20 seconds between the peaks 
in buoyancy term {q^^^ ) and turbulent kinetic 
energy. This is about half the time it takes 
the flow to transit a distance icz, the depth 
of the convection zone. It also corresponds to 
an e-folding time for turbulent kinetic energy 
decay due to Kolmogorov damping, where 

Ek = U^ms/^d = {ul^j2){2Urms/U)- Power 

spectra for both variables peak at 89 seconds; 
an average transit time is 51 s. 

Figure [1] shows multimodal behavior in 
the 3D simulations; prelim inary results from 



a qu antitative analysis (jMeakin &: Arnett 



201 ll ) using principal orthogonal decomposi- 
tion (POD) indicates that a single dominant 
mode has about 43% of the kinetic energy, the 
first five modes have 75%, and 90% is reached 
with the eleventh mode. There is a strong 
dominant mode but also significant energy in 
several other modes; the modes interact in a 
nonlinear and dynamic way. 

The buoyancy power is a large scale fea- 
ture and is strongly anisotropic (plumes move 



largest eddies. A more careful discussion of the veloc- 
ities in terms of principal orthogonal decomposition is 
in preparation. 



vertically). The dissipation implied by the 
turbulence occurs at the Kolmogorov scale 
(which is tiny); this dissipation is wide-sp read 
in space ( Arnett. Meakin. &: Young 20091 ). in- 
cluding the entire turbulent region on average, 
and bounded by the stably stratified layers. 
Because buoyancy and dissipation occur at 
vastly different length scales, they are weakly 
coupled. 

In low Mach number flow, on average over 
time, the buoyant driving must balance the 
turbulent d amping for a quasi-steady state 
to exist (see lArnett. Meakin. fc Young! (|2009l ) 
for a detailed discussion, especially their 
Eq. 33). The time averaged viscous dissi- 
pation is 



^turb — (lintl^CZ 



^rms 



'/£d-lMu^/£d, (9) 



where id is the dissipation length, icz is 
the depth of the convection zone, and the 
over-lines indicate an average over time (two 
turnover times or four transit times, about 
200 seconds, in the analysis). 

From Equation [9] it is clear that the turbu- 
lent dissipation is of third order in the tur- 
bulent velocity, while the buoyancy power, 
q = —gu'p'/po, is second order (because p' 
is first order). This means that in the turbu- 
lent regime there is a unique turbulent veloc- 
ity which satisfies the condition that buoyancy 
power balance turbulent damping. This has 
the nature of an eigenvalue problem. The as- 
sumption of a turbulent cascade implies finite 
turbulent velocities, with macroscopic struc- 
ture. Equation [9] applies as an integral over 
a region much larger than the Kolmogorov 
scale; at the Kolmogorov scale microscopic 
dissipation occurs by the Navier-Stokes term 
dLandau fc Lifshitd (Il959l ). §15). If we dot 
this Navier-Stokes term with the velocity to 
generate a rate of kinetic energy destruction 
by viscosity (see Equation [5]) , we have 



{u-du/dt)i 



z/u-V^u 



(10) 
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which does scale as the velocity to second or- 
der. When integrated over the turbulent cas- 
cade, Equation [10] becomes Equation [9l The 
extra factor of velocity comes from the rate at 
which turbulent kinetic energy is delivered to 
the Kolmogorov scale. Turbulence is an inher- 
ently nonlinear process, which sets up its own 
scales and structures. Care must be taken in 
analyzing it with techniques involving expan- 
sion in series [^. 

The dissipation is seen to be driven by the 
isotropic part of the flow ut, which accounts 

for about three-fourths of the kinetic energy 

accord ing to estimates bv lArnett. Meakin. Sz Young 
(|2009l ) (which are roughly similar to those 



from POD analysis by iMeakin &: Arnett 
( 201ll )). Note that v'^^^ = which gives 
the vertical offset between the solid and dot- 
ted curves in Figure [TJ Consider the average 
conditions from 200 to 800 seconds in the 
simulations. The average level of buoyancy 
(actually q^^^^ is plotted) is higher than that 
of kinetic energy in isotropic turbulence. If we 
equate (1) the power generated by buoyancy, 
to (2) the dissipation due to turbulence, both 
averaged over a few transit times, we find 



cz 



int 



0.85. 



(11) 



Thus, in the iMeakin &: Arnett (l2007bl ) simu- 
lations, the dissipation length is found to be 
essentially the depth of the convective zone, 
consistent with Kolmogorov theory. 

These results are far more general than the 
particular stellar situation we have discussed. 
While the neutrino cooling may seem exotic to 
some stellar evolutionists, in fact it behaves 
somewhat like the more familiar cooling by 
radiative diffusion, and has no strange effect 
on the turb ulence. We note that the original 
simulations ( Meakin &: Arnett 2007bl ) . which 



''See I Arnett fc Meakinl (|2011bt ) for an example of the 
danger of using linear stability analysis for turbulence 
in stars. 



included core hydrogen burning cooled by ra- 
diative diffusion as well as oxygen shell burn- 
ing cooled by neutrinos, explicitly showed this 
similarity. 

In the long term, the thermal state of the 
convection zone is supposed to evolve toward 
a global thermal balance between total heat- 
ing by nuclear bur ning and total coo ling by 
neutrino emission ( ArnettI 1972 . 19961 ). This 
is illustrated in Figure [21 which shows time 
average luminosities of heating and of cooling, 
as a function of the entropy of the convection 
zone. At lower entropy, as in the 3D simu- 
lation, the heating is dominant, causing the 
entropy to increase. This is accomplished pri- 
marily by expansion, with a small increase in 
temperature. The decrease in density causes 
a bigger change in nuclear burning than in 
neutrino cooling, so that a thermal balance 
would be attained (shown as a vertical line) 
when the heating and cooling curves cross. 

The gradual rise in turbulent kinetic en- 
ergy in the simulation, shown in Figure [H oc- 
curs in a shell which is below the entropy of 
balanced heating and cooling, so that heat- 
ing dominates. The pulses are much faster 
than this secular evolution, which changes on 
a time scale of t > 10^ seconds. 

4. The Lorenz Solution 

The Lorenz model is a convective roll, or 
cell, whose dynamics are described by three 
amplitude equations. This is a simple ex- 
ample of an elegant method of reduction of 
turbulent flow to a low-order set of dynam- 
ical equations using amplitudes of a proper 
orthogonal decomposition (POD) of numeri- 
cal si mulations or extensive experi mental data 
sets ( Holmes. Lumlev. Sz Berkood 1996 ). 

The Lorenz model is a better mathemati- 
cal representation of the dynamics of a con- 
vective cell than MLT in that the accelera- 
tion and deceleration over a convective cy- 
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cle are integrated to determine the motion, 
rather than prescribed. Because it uses a sin- 
gle mode, the Lorenz model does not have sen- 
sitivity to local variation; it is a global model, 
unlike MLT, which is local. Lorenz devised 
the model as a point of principle test of me- 
teorological convection, which like the stel- 
lar problem, is damped by a turbulent cas- 
cade. To make us e of this extensive literature 
dCvitanovid Il989l: iF risch' 'l995|; Icieickl Il987l : 
Thompson &: Stewar d 1986), we use the orig- 
inal version of iLoren z (19631), with the same 
Prandtl and Rayleigh numbers, to explore the 
implications on dynamics of the strange at- 
tr actor. The original Lorenz formulation is 
a 2D, low Mach number, and single mode 
model (equivalent to Figure [3]). A transi- 
tio n to 3D may be made using the solutions 
of Chandrasekhail ( 196 ll ). but is not neces- 
sary for present purposes . Real convection 



(Libchaber Mauer 



19821 ) is expected to be 



single mode only near the onset of convective 
instability. 

The Lorenz equations (see Appendix A) 
are: 




CZ Entropy 



Fig. 2. — Nuclear heating and neutrino cool- 
ing in an oxygen burning shell, as a function 
of entropy. At lower entropy, as in the sim- 
ulations, nuclear heating dominates because 
of its density dependence. For the same tem- 
peratures but higher entropy, the inhibition 
of neutrino-antineutrino emission is reduced. 
The shell entropy slowly increases until cool- 
ing balances heating, at which point thermal 
balance is attained. 



dX/dr = -aX + aY (12) 
dY/dr = -XZ + rX-Y (13) 
dZ/dr = XY - bZ, (14) 

In the classical Lorenz model, X is propor- 
tional to the speed of convective motion, Y is 
proportional to the temperature difference be- 
tween ascending and descending flow, and Z 
is proportional to the distortion of the vertical 
temp erature profile from adiabatic (jLorenz 
19631 ). Here r is a time in thermal diffusion 
units, a is the effective Prandtl number, r is 
the ratio of the Rayleigh number to its critical 
value for onset of convection, and a is the ra- 
tio of the depth to the width o f the convectiv e 
cell, so that b = 4/(1 + a^); see Lorenz ( 19631 ). 



The Prandtl number is the ratio of coefficients 
of the viscous dissipation term to the thermal 




Fig. 3. — The Lorenz Model of Convection: 
Convection in a Loop. 
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mixing term. LorenqJ took a = 10, and chose 
the most unstable mode so that = 1/2, and 
b = 8/3. For this mode, and a Rayleigh num- 
ber of r = 470/19 = 24.74 times the critical 
value, steady flow becomes unstable. 

Figure [3] illustrates the structure and nota- 
tion in the Lorenz model. The Prandtl num- 
ber may be expressed as a = T^ad/^visc, the 
ratio of the radiative cooling time scale to the 
viscous time scale. The viscous damping time 
is taken to be constant, Tyisc — 

1/r. The 

Rayleigh number, in units of its value at the 
onset of laminar convection, is 

r = g^TTi/iTKTo, (15) 

where K = 1/Trad, f^T = -{dhip/dhiT)p, 
and g is the gravitational acceleration. The 
lapse rate of the adiabatic background Te is 
g/Cp, so Ti = g{l/2)/Cp. In astrophysical 
notation, 

V - V, = Va{T2/Ti - 1), (16) 

where T2 is the temperature variation in the 
vertical direction (Fig. [3|). The depth of the 
roll, in pressure scale heights, is 

i/Hp = ln[(ro + Ti)/(ro - Ti)]. (17) 

V a 

4.1. Energy fluxes 

In the meteorological case the fluid motion 
(wind) is of primary interest, however the en- 
ergy fluxe^ provide a means to connect with 
stellar observations. 

justification is simply that this allows the three 
Lorenz equations to capture chaotic behavior. See Ap- 
pendix B for a discussion of the physical implications 
of the Prandtl number for stars. 
*In discussion of the Lorenz model, both heat flux and 
enthalpy flux are mentioned; which is correct? De- 
pending upon the physics perspective, both may be. 
In the simplest Lorenz model, the fluid in strictly in- 
compressible, so that the work term PdV is strictly 
zero, and in this case only heat content can carry en- 
ergy, so heat flux is relevant. However, we may inter- 



The vertical enthalpy flux due to radiative 
diffusion of the background is 

Fe = -pvTCpdTE/dz (18) 
= puTCpTi{2/t), (19) 

which is constant in the Lorenz model. 
The thermal conduction coefficient is vt = 
{AacT^/pK)/{pCp) and K = i/t(1/^)^- Here 
z = — I cos(/), the radius of the roll being £/2. 
Because the background temperature is taken 
to be linear in z, the divergence of this flux is 
identically zero, and therefore gives no local 
heating or cooling. The additional vertical 
enthalpy flux, due to the temperature pertur- 
bation, is 

F, = -puTCpd(T-TE)/dz (20) 
= piyTCpiT2-Ti)i2/e). (21) 

If we define a potential temperatur^ T4 = 
T1—T2, the vertical flux is separated into two 
parts, so that the adiabatic background value 
is denoted by Ti and the changes caused by 
motion are contained in T4. 

The flux in the horizontal direction is 

Fy = -pvTCpdT/dy (22) 
= -pvTCpT:,{2/i), (23) 

which averages to zero by symmetry; here y = 
I sin (/>. Both F^ and Fy are proportional to a 
potential temperature which varies in time. 

The net vertical enthalpy flux will not be 
constant in height z, so that its divergence is 
nonzero. This implies local heating/cooling. 



pret the equations in terms of a stratified system with 
low-Mach number flow, having a stratiflcation in den- 
sity as well as temperature. Then dV is not zero, PdV 
work is done, and the relevant flux is the enthalpy flux. 
Here potential tempe rature is used as in fluid dy nam- 
ics and meteorology (|Duttonlll98^ : lTrittonlll988l l: the 
potential temperature is measured relative to an adi- 
abatic background, and may be small even if the con- 
vective region is strongly stratifled (as is often the case 
in stars). 
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which would have to be compensated for in a 
steady state, and is a consequence of consid- 
ering only a single mode. Smaller scale modes 
would be neede d to deal with the loc a l heat - 
ing/cooling; see ICanuto &: MazziteUi 
Canuto et al. I (j 19961 ) for efforts to include a 
full spectrum of modes. To get the average 
flux through the cell for a single mode, we 
take the double projection, of u and of T on 
the vertical direction. 



pUzCpT 
pCpuT^ SIT? 



(24) 



The enthalpy flux is zero at the top and bot- 
tom of the roll [(f) = tt and 0), and a maxi- 
mum at the midpoint {(f) = 7r/2). The average 
value of sin'^ (f) over the roll, ^ = to 27r, is 
0.5. This is to be compa. r ed wi th the q^; of 
Arnett. Meakin. k Youn"g| (l2009l l. which was 
the range 0.68 to 0.85 for the simulations then 
available. Because the simulations are multi- 
mode, this slightly higher value for the corre- 
lation seems natural. 

Integrating over the cell, < (/> < 27r, gives 



= \pCpuT?^ 



(25) 



which is the vertically averaged enthalpy flux 
of the cell. This is the "convective flux" in 
stellar models. The up-flow enthalpy flux 
equals the down-flow value, and both are pos- 
itive (hot up-flows and cold down-flows both 
transport energy upward). 

The ratio of enthalpy flux to the total ra- 
diative flux {Fz + Fe) is 



FjFr = {1/Aut){uT^/T2) 



(26) 



which is a function of time. 

The kinetic flux exactly cancels in this for- 
mulation, as in MLT, with the up-flow being 
the negative of the down-flow. This is not 
true in the simulations. The kinetic flux does 
approach zero for convective regions which 



are so thin that they are almost unstrati- 
fied. However, stratified convection zones 
have an asymmetry in up and down flows, giv- 
ing a modest net kinetic energy flux (down- 
ward for pure top driving, up for pure bot- 
tom driving, and both for more general cases, 
Meakin fc ArnettI (|^09|)). 

We may express fluxes in units of the radia- 
tive flux of the background, Fe, and in terms 
of the variables in the Lorenz equations. The 
excess vertical radiative flux is 



Fz/Fe = n/Ti 
= Z/r, 



(27) 



where Ti = g£/2Cp, Tq = gHp/{Cp - Cy) 
and 7 = Cp/Cy, and 



7-1/ gQ 



27 \HpKT 
The horizontal radiative flux is 

Fy/FE = n/Ti 
= Y/r. 



(28) 



(29) 



The enthalpy flux of the cell, in the vertical 
direction, is 



Fe/FE = XY/2r. 



(30) 



In general, without need for any additional 
mechanism to cause variability, the net flux 
of energy through a turbulent cell varies with 
time. 

4.2. Steady-state solutions 

A steady state solution for the Lorenz equa- 
tions is 

X = y = Z = 0, (31) 

and, for r > 1, a second solution, for a stable 
convective roll, appears: 

X = Y = ±[h{r -l)Y''^;Z = r -I. (32) 
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The second solution is unstable if o" > 2 and 
r > rc where 



a{a + b + 3)/{a-b-l) 



(33) 



and this in stability give s rise to fluctuations 
in velocity ( Lorenz 19631 ) . The instability also 
gives fluctuations in energy flow through the 
cell. 

The steady state solutions might be of 
practical value for stellar evolutionary codes, 
to the extent that they provide an estimate of 
the average behavior of the convective cell. 
However, the fluctuations are large (in no 
sense are they "small perturbations"), and 
may cause nonlinear complications. For ex- 
ample, a thermonuclear runaway would be 
sensitive to the largest value of the tem- 
perature fluctuations bec ause that would 
affect the burning rate (jArnett &: Meakin 
201 Ibl ^. The fluctuations drive entrainment 
episodes, and affect the ni ixing of composition 
(jMeakin k ArnettI lioOTbl ) . The fluctuations 
may be able to modify the driving of pulsa- 
tions in stars with vigorous convection zones. 
In red supergiants, coupling of turbulent fluc- 
tuations in the surface convection zones to 
both radial and non-radial pulsations seems 
likely (see below). 

4.3. Non-steady Solutions 

In order to compare the Lorenz model to 
3D simulations, it is desirable to have compa- 
rable starting conditions. Arbitrary choices 
can give large initial transients before the at- 
tractor controls the behavior. The 3D simu- 
lations start with zero velocity but finite tem- 
perature deviations from an adiabatic gradi- 
ent. We take 



X = 0,Y = [b{r - 1)]2,Z 



1, 



(34) 



as initial conditions, which sets the velocity 
to zero and uses steady state values for the 
temperature fluctuations. 



Figure [T] shows the behavior of the Lorenz 
model of convection (panel labeled LORENZ), 
for a similar number of pulses as shown for 
the 3D simulation, and the same variables: 
buoyancy power XY/(b{r — 1)) to the two- 
thirds power, and kinetic energy per unit mass 
X"^ /{b{r — 1)). The factors b{r — 1) are cho- 
sen so that the steady state values are nor- 
malized to unity, but are qualitatively cor- 
rect for comparison with the 3D simulations. 
The velocities in Fig. [T] are not precisely the 
same, the 3D simulations giving turbulent ve- 
locity while the Lorenz speeds are more ap- 
propriate to "coherent structures". They are 
related (see § 3.1), as the large scale ve- 
locities become unstable and turbulent, and 
their kinetic energies are the same order of 
magnitude (Arnett. Meakin. fc Young 20091 : 



Meakin fc ArnettI bo 111 ). 



The time is measured in the dimension- 
less units of Lorenz; a turnover time is two 
transit times, and roughly the time between 
peaks. We see that the peaks in buoyancy 
power slightly precede those in kinetic energy, 
as in the 3D simulations but less dramatically. 
This difference is related to the fact that the 
Lorenz model has viscous dissipation acting 
directly on the large scale velocities, while the 
3D simulations have dissipation at the Kol- 
mogorov scale, which is separated from the 
large scale by the turbulent cascade, involv- 
ing many, many reductions in length scales 
(see Section IB.3|) . Additional modes would 
fill in the "valleys" in the Lorenz model (see 
Over the time shown, the average kinetic 
energy is 0.968 of the formal steady state solu- 
tion, which is shown in Figure [1] as the dashed 
horizontal line in the Lorenz panel. 

The degrees of freedom in the three- 
dimensional simulations, and in the Lorenz 
model, are dramatically different. The floating- 
point operation count differs by a large factor: 
~ 10® (several times 10® zones times 7 effec- 
tive scalar variables for simulations, versus 
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three amplitudes for the Lorenz model) . With 
such an extreme compaction, it is striking that 
they give a similar picture for fluctuations in 
convection. 

4.4. Long Term Behavior 

Unlike the full 3D simulations, it is triv- 
ial to extend the Lorenz model to later times. 
Figure H] shows an extension in time by a fac- 
tor of 10. There is a relatively steady growth 
in amplitude up to time ~7, at which point 
chaotic behavior begins to appear. The fluc- 

2 

tuations in kinetic energy (and q?^^) are large, 
sometimes exceeding the steady state value 
(~1) by a factor of 4. The convective lumi- 
nosity from a single Lorenz model shows the 
same qualitative behavior. 

This drastic behavior raises two interesting 
possibilities: 

1. The numerical simulations will follow 
the solutions of the Lorenz equation, 
and exhibit vigorous and chaotic fluc- 
tuations at later times, or, 

2. The multimode behavior of the simula- 
tions will allow even stronger dissipa- 
tion, preventing extreme behavior. 

Either way, the result is important for the evo- 
lution of supernova progenitors, especially re- 
garding the effects of fluctuations on mixing 
(yields) and outbursts. Simulations of mul- 
tiple burning shells (C, Ne, O and Si) in 2D 
have been continued further than the 3D sim- 
ulations for the oxygen burning shell. They 
appear to follow the flrst option, so that the 
prediction from 2D simulations is that core 
collapse progenitors will h ave violent erup- 
tions p rior to core collapse (jArnett &: Meakin 
2011br). Pull 3D simulations need to be ex- 



4.4- 1- Duration 

The oxygen shell might not last for the 
~70 transit times shown in Figure HJ a lin- 
ear estimate suggests that it consumes its fuel 
in about 100 travers al times (the Damkohler 
numb er is Da < O.Ol. lArnett. Meakin. &: Young 
(Hooi)), so the background evolution should 
not be neglected for such time intervals. 
Oxygen burning is likely to be a more dy- 
na mic event than previously s uppo sed (e.g., 
by Wooslev. Heger. Weaver ( 20021 )). Aver- 
aging over multiple cells may moderate the 
net fluctuation (see below). 

For convection zones in other stars, the 
number of traversal times available may be 
much larger. The sun has a deep surface 
convection zone (20 pressure scale heights); a 
plume would, if unimpeded, fall through the 
convection zone in about 2 hours; it would 
take about 7 centuries to process all the mass 
of the convection zone through the surface lay- 
ers. These times bracket the effective mixing 



tended to later times at which the Lorenz 
model predicts chaotic behavior. 



time, so that of order 10 to 10 mixing times 
would occur over the age of the present-day 
sun. 



4-4-^- Chaotic Behavior 

Figures [5] and [6] show the fam iliar long term 
behavior of the Lorenz model (jGleickl 119871 ) . 
Figure [5] shows the trajectory in X and Y 
space for the flrst 70 dimensionless time in- 
tervals. X is the dimensionless speed and 
Y the dimensionless temperature variation in 
the horizontal direction. Both X and Y switch 
signs, but there is an overall correlation be- 
cause there is a net enthalpy flux, which is 
proportional to XY. The correlation is in the 
sense that hot regions rise and cool regions 
sink. The radiative flux associated with the 
horizontal temperature variation is propor- 
tional to Y, and while on average over long 
times is zero, it achieves this in segments of 
time in which first one then the other direc- 
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Fig. 4. — Later Behavior of the Lorenz Model 
of Convection, for a = 10, r = 28, and 
b = 8/3. The format is the same as Figure [H 
but the time interval is longer by a factor of 
10. Chaotic behavior is beginning to appear, 
with fluctuations in kinetic energy and buoy- 
ant driving sometimes exceeding the steady 
state value by a factor of 4. 




Fig. 5. — Long Term Behavior of the classic 
Lorenz Model for X (speed) and Y (horizontal 
temperature fluctuation), for a = 10, r = 28, 
and b = 8/3. The trajectory is shown for 70 
Lorenz time units. Both X and Y switch signs, 
but there is an overall correlation (hot zones 
rise while cool zones sink). This correlation 
gives rise to positive values for buoyancy flux 
and enthalpy flux. 



tion (sign) are favored. This is the strange 
attractor at work. 

Figure [S] shows the same time interval, but 
plots Y and Z (the dimensionless temperature 
in the vertical direction, parallel to the grav- 
ity vector). Notice that Z has positive val- 
ues about which it orbits. The stratification 
breaks the symmetry in vertical versus hor- 
izontal directions. The radiative flux in the 
vertical direction has a steady part and a fluc- 
tuating part. The latter is proportional to Z, 
and is positive (outward energy flow). 

In MLT, the Schwarzschild discriminant 
{S = V — Va) is the temperature excess above 
adiabatic; this excess is assumed to be pro- 
portional to the enthalpy flux. This is incon- 
sistent with the Lorenz model, because the 
velocity X can have both signs while Z does 
not. The error arises because in MLT the 
speed of convection is taken to be intrinsically 
positive (to avoid this problem), and may be 
traced back to a lack of conservation of mass 
("blobs dissolve" into the environment rather 
than flow back to complete a cycle). 

4.4.3. MLT 

What happens if we reduce the three equa- 
tions of Lorenz to two, forcing one variable 
to be at its steady state value? This is simi- 
lar to the MLT approach (assuming the mix- 
ing length parameter is chosen so that the 
kinetic energy scale is physically correct; see 
Arnett. Meakin. &: Young] ( 20ld )). Enforcing 
the steady state value for the vertical tem- 
perature excess, Z = r — 1, but allowing 
X and Y to vary, corresponds to a model 
with a single temperature variable (like MLT). 
Such integrations are shown in Fig. [71 There 
are no pulses in kinetic energy or buoyancy; 
the curves quickly approach a constant value. 
Convection proceeds by steady motion in a 
roll; a finite XY is required to give torque to 
make the roll. This "two equation model" no 
longer has a strange attractor; the pulses have 
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Fig. 6. — Long Term Behavior of the clas- 
sic Lorenz Model, for a = 10, r = 28, and 
b = 8/3. Same as Figure [5l but with Y and 
Z (vertical temperature fluctuation) shown. 
Z is a positive for these parameters. Y has 
both positive and negative steady state val- 
ues, about which it orbits. 




Fig. 7. — Initial Behavior of the two Equation 
Model of Convection, for a = 10, r = 28, and 
b = 8/3. The critical value for instability of 
the convective roll is Tc = 24.74. The format 
is the same as Figure [U The turbulent fluc- 
tuations are eliminated, as they are in MLT. 



been eliminated. This explains why stellar 
evolutionary calculations which use MLT do 
not show these fluctuation^^. 

5. Cells and Shells 

5.1. Multiple Modes in Cascade 

Figure [1] shows that the primary differ- 
ence in the Lorenz model and the 3D sim- 
ulations is that the Lorenz model has only 
a single mode, while the simulations are ob- 
viously multimodal. This difference may be 
superficial. The Lorenz model in this appli- 
cation (as in the original meteorological one) 
has additional modes implied by the turbu- 
lent cascade which mediates the damping (i.e., 
they are implicitly in F; see Appendices ^A] 
and ^B]). a simple Richardson cascade was 
discussed in §B.3( in which /, the fractional 
change in length scale for each ste p in the cas- 
cade, is assume d to be a constant ( FrischlE995l : 
Davidsonll2004l ) . This is not very plausible for 
the largest scale modes because they are the 
most sensitive to boundary conditions (they 
must fit into the convective region) , but is sim- 
ple and instructive. The fractional time spent 
in the cascade for each mode may be shown to 
be the fractional kinetic energy in that mode. 
Using ffHH this gives /(2/3)(n-i) - /2/3) for 
n = 1,2,---, or roughly 0.37, 0.23, and 0.14 
for the first three. There is a dominant mode 
accompanied by several weaker but significant 
ones. 

One way to proceed would be to introduce 
additional modes into the Lorenz model (cho- 
sen with guidance by the 3D simulations), 
and to generate an expanded set of amplitude 
equations which generalize the three of Lorenz 
(Eq.iniini andlM]). This would make the sys- 
tem multimodal, and allow for modal interac- 



^" In analogy, the two-body (Kepler) problem in celes- 
tial mechanics is well behaved, while the unres tricted 
three -body problem is far more complex ((Poincarel 
Il893l l. 
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Toroidal Flow 




Up-draft 

Fig. 8. — The Lorenz Model extended: con- 
vection in a sphere composed of a ceh of 
toroidal shape. With no rotation to break the 
symmetry, the direction of the upflow vector 
is not constrained, and will be chaotically cho- 
sen. 
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Fig. 9. — The Lorenz Model extended: Con- 
vection in a shell composed of cells. Notice 
the alternation of the sign of rotation. This 
may be thought of as a cross sectional view of 
infinitely long cylindrical rolls, or of a set of 
toroidal c ells, with pairwise alt ernating vor- 
ticity (see IChandrasekhaii (|l96ll ). §16). Each 
cell may exhibit independent fluctuations in 
time and space. 



tions. This approach is a simplified version of 
an el egant proposal by Lumley and colla bora- 
tors ( Holmes. Lumlev. &: Berkooz 19961 ): em- 
pirical eigenfunctions are constructed by POD 
from simulations, introduced into the differ- 
ential equations to derive a set of nonlinear 
ODE's for the amplitudes of each mode (a 
Galerkin projection), and this set of ODE's is 
solved to generate the evolution of the aver- 
age prope rties of the turbulence. This is being 
explored (jMeakin &: Arnettll201ll ) as a way to 
capture more fully the dynamics and multi- 
modal behavior seen in the 3D simulations. 

5.2. Multiple Cells 

Suppose we envisage the convective re- 
gion to be populated by cells, each of which 
is a separate Lorenz model representing the 
largest mode, a convective roll, in that cell. 
A special case is the convective core: the ge- 
ometry is not that of a layer, but a spherqlj, 
as is indicated in Figure [8l The flow has a 
toroidal structure. In this case the gravita- 
tional acceleration goes to zero at the origin 
(the center of mass of the star), but the veloc- 
ity need not be zero there. If there is no net 
rotation, then the direction of the up-flow is 
undetermined, and will be chosen chaotically 
by the turbulent flow. 

For a convective shell, we imagine that the 
layer is filled by cells, as shown schematically 
in Figure El For the oxygen burning shell, 
the inner and outer radii are about 4 x 10® 
and 8 x 10® cm, respectively. The area of 
the spherical shell, evaluated at the midpoint 
in radius, is 47r(6 x 10®)^, and the cell is 
taken to be roughly square, so that its area 
is i"^ = [A X 10®)^. The ratio of sheh area 
to cell area is about 9tt, so we assume there 
are roughly 30 cells spread over the spheri- 



^ IWoodward. Porter, fc Jacobs! (|2003l ) found such a " si- 
ant dipole" behavior in their 3D simulations of almost 
fully convective spheres. 
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cal shell. In general, the number of cells in 
a shell will depend upon the geometry of the 
convection zone. 

If the cells are not synchronou^^, but act 
independently, the effect of the pulses will 
be smoothed when summed over the whole 
shell. However, the cells may interact con- 
structively; the solution to this more complex 
problem remains opeij^. At issue are both 
the interactions between cells in a single con- 
vective region, and interactions between mul- 
tiple convective regions associated with differ- 
ent burning shells. 

The individual cells exhibit fluctuations 
not only in time, but also in space. Each cell 
represents a mode which is unstable, and de- 
stroys itself and reforms, usually somewhere 
else. Because the medium is fluid, the pat- 
tern of cells is much more dynamic and less 
regular than that of a crystalline solid, so the 
Fig. should be interpreted as representing a 
snapshot of a system which fluctuates in both 
space and time. Hinode/SQT observa tions of 
the solar surface ( Berger. et al. Il20ld l reveal 
such a highly dynamic, complex structurJ^. 

5.3. Irregular Variables 

Figure [TU] explores an idea o f Martin 
Schwarzschild ( Schwarzschild 19751 ). who es- 
timated the number of convective cells in the 
sun and in red giants and supergiants. He ar- 



^^See lAlligood. Sauer fc Yorkd l| 19961 ) for a discussion of 
synchronization of chaotic orbits. 

" Figure 23 in iMeakin fc ArnettI (|2007bl ) suggests the 
complexity of the cell interaction within a single con- 
vective region. The original simulations were on a 
wedge, of 27° in theta and in phi. Simulations with 
larger aspect ratio (larger angle wedges) do show a 
moderation of the total fluctuation, in qualitative 
agreement with the dis cussion above. Figure 2 in 
iMeakin fc Arnet3 (|2006l ) indicates the complexity of 
interaction between multiple burning regions even in 
2D. 

^''We thank Dr. A. Title for providing a copy of his 
spectacular movies of data and simulations. 



gued that only a modest number of cells (~12) 
would exist at any given time in a red giant or 
supergiant. To illustrate the point, we use the 
Lorenz model to approximate the behavior of 
a convective cell. We have estimated the con- 
vective flux for 12 cells at random phase, by 
adding 12 time sequences from an computa- 
tion like Figure S] (but extended to t=800), 
starting at 12 randomly chosen times in this 
interval. This time sampling is intended to ap- 
proximate a spatial ensemble average. Their 
flux, summed and normalized, is shown 
solid line, for a dimensionless time from 60 to 
70, which corresponds to roughly 20 pulses. 
The signal is noisy and looks "chaotic". 

These fluctuating cells make up a convec- 
tive region, and will couple to the normal 
modes of the star to cause both radial and 
non-radial pulsations. The amplitude of these 
pulsations will depend upon the overlap inte- 
grals between the normal modes and the cell 
motion, and the stellar damping. This sug- 
gests that the noisy behavior will be combined 
with the relatively cleaner periodicity of the 
normal modes, giving a power spectrum with 
a base like that shown in Figure [TTl but with 
superimposed spikes corresponding to the ex- 
cited normal modes. While turbulent convec- 
tion alone is sufficient to cause luminosity fluc- 
tuations, it occurs in regions of high opacity 
and partial ionization, which also drive pulsa- 
tion, so that composite behavior and multiple 
periods may be expected. 

Joel Stebbins (pioneer of photoelectric as- 
tronomy) monitored the brightness of Betel- 
geuse (a Orionis) from 1917 to 1931, and 
concluded that "there is no law or order in 



the rapid changes of Betelgeuse" ([Goldberg 



1984 ). which seems apt for Fig. [10] (the Lorenz 
strange attra ctor) as well. More modern 



observations ( Kiss. Szabo. Bedding! 20061 ) 
show a strong broad-band noise component 
in the photometric variability. The irregular 
fluctuations of the light curve are aperiodic. 
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and resemble a series of outbursts. Direct 3D 



Fig. 10. — Fluctuations of Luminosity in Con- 
vective Layer of 12 Cells of Random Phase, 
for a = W, r = 28, and b = 8/3. The 
dimensionless flux (luminosity) is shown for 
a convective layer with 12 visible Lorenz 
cells. The luminosity variations are large 
and seemingly chaotic, suggestive of irreg- 
ular variables and Betelgeuse in particular 
(|Kiss. Szabo. fc Beddin jl2006l l. 




Fig. 11. — Power spectrum of the luminos- 
ity fluctuations implied by turbulence alone. 
There is no sharp peak, but a broad distri- 
bution of power, as would be expected from 
a chaotic source. Resonant interaction with 
normal modes of the star could add peaks to 
the spectrum, which would provide an obser- 
vational constraint on interior structure, anal- 
ogous to astro-seismology of solar-like stars. 



simulations of Betelgeuse (IChiavassa. et al 



20101) show the same complex behavior (and 



have the advantage that they predict the de- 
tailed spectral behavior as well). This should 
be no surprise; the 3D equations have embed- 
ded in them a strange attractor. 

Increasing the number of cells reduces the 
level of fluctuation about the mean. Averag- 
ing over the 2 million granules of the sun gives 
a very stable luminosity, which would plot as 
a straight line in Fig. [10] (even 2000 random 
cells do this). However, the size of the cells at 
the bottom of the solar convection zone will be 
larger (hence fewer cells) , and if chaotic might 
give a long term modulation to the solar lumi- 
nosity. Full star simulations of the whole so- 
lar convective zone, with sufficient numerical 
resolution to give well developed turbulence, 
should shed light on this issue. 

Application to turbulent stellar atmo- 
spheres requires surmounting two difficulties: 
(1) the flow is no longer low-Mach number 
(see § 2.2), and (2) the ionization zone causes 
dramatic changes in opacity (assumed con- 
stant in the Lorenz model). Fortunately 3D 
atmospheres exist, and analysis such as we 
have done on stellar interior convection simu- 
lations is feasible. 

While this paper was in preparation, 
Stotherd (j2O10l ) re-examined the idea that gi- 



ant convective cell turnover is the explana- 
tion of the long secondary pe riod observed in 
semiregular red variable stars (jStothers &: Leung 
197ll ). including Betelgeuse and Antares. 



Stothers used MLT to derive a velocity scale 
for the overturn, also r elying on general fea - 
ture of simulatio n s of Chan fc Sofia ( 19961 ): 



see 



3 in IStothersI ([20101). This theory ap- 



pears to work directly clS 3) complement to 
the discussion above. The use of a Lorenz 



^ lArnett. Meakin. &: Yound (|2009l ) discuss and contrast 
these and other simulations. 
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model for the giant cells already implies real 
dynamics. The strange attractor necessarily 
provides vari a.bility in luminosity, with its own 
quasi-period ( Wood. Olivier &: Kawalei 20041 ') 
and velocity scale (see §4 in IStothersI (|201ol V). 
Several giant cellq^ are at work, each with a 
quasi-period of order of the transit time, and 
therefore similar to the estimate of Stothers, 
and the observations. The larger convective 
velocities needed are simply a necessary con- 
sequence of the dy namics implied by the con- 
vecti ve luminosity (jArnett. Meakin. &: Young 
2OIOI ). 



The introduction of turbulence as an active 
agen t in the discus si on of stellar variability 



Stothersl(|2010l 'l: IWood. Olivier fc Kawaler 



(|2004l )) seems timely. An interesting improve- 



ment would be to use POD ei npirical eigenval- 



ues fr om simulations (e.g., IChiavassa. et al, 
( 20ld ). § 15.11 above), and develop a low or- 
der dynamical model to explore long duration 
behavior. 

6. Summary 

We have identified a major new feature of 
stellar physics: chaotic behavior due to tur- 
bulent fluctuations in stellar convection, and 
corresponding luminosity fluctuations. While 
the simulations upon which the analysis was 
based were fully compressible, the theory uses 
the approximation of sub-sonic flow. Both nu- 
merically and analytically, a strange attrac- 
tor like that of Lorenz ( 19631 ) seems to appear 
naturally in stellar convection. 

As a first approximation to more rigorous 
analysis, we have applied the Lorenz model 
to kinetic energy fluctuations in the oxygen 
burning shell, to the turbulent energy cascade, 
and to fluctuations in luminosity in irregular 



This idea can be tested observationally and numeri- 
cally; the structure of the convection region will con- 
strain the number of cells, which may be compared 
with the amplitude of luminosity fluctuations. 



variables. 

Figure [1] shows a comparison of the tur- 
bulent kinetic energy fluctuations in 3D sim- 
ulations of turbulent flow and in the Lorenz 
model. No parameters were adjusted to give a 
fit. Additional modes, appropriate for turbu- 
lent flow, would improve the comparison fur- 
ther. 

This suggests a new, inherently nonlin- 
ear mechanism for variability in stars, the 
r-mechanism, which is caused by luminos- 
ity fluctuations directly associated with tur- 
bulent convective cells. Because the mech- 
anism is nonlinear, it is not captured by 
linear stability analysis, which is a main- 
stay of variable star theory. Such lumi- 
nosity fluctuations may have been observed 
already in the broad-band noise se en in 



Betelgeuse (alpha Orionis (iGoldbera Il984l : 



Kiss. Szabo. &: Bedding] |2006| ). and in the 



long secondary periods in pulsating AGB stars 
(|Wood. Olivier Kawalel2004l : IStotherdl20ld l , 
and are expected to be observable in principle 
in all stars with extensive surface convection 
zones, including those with "solar-like" vari- 
ability. This mechanism is probably the cause 
of the strongly driven pulsa t ions found by 
Woodward. Porter. &: Jacobs! (I2OO3I I in their 



"red giant" model; the development of those 
large pulsations was a clue which may now be 
more fully understood. 

Such fluctuations provide a source of per- 
turbations for instabilities, and may induce 
mixing not presently accounted for in stellar 
evolutionary calculations. The fluctuations in 
convective velocity are comparable to average 
values. If these fluctuations couple to nuclear 
burning, as for example in cases of degenerate 
ignition, shell flashes, or later stages of oxygen 
and silicon burning^, outbursts may develop. 



^^For new developments since th is paper was submitted, 
see lArnett fc MeakinI (j2011bh . where it is suereested 
that core collapse progenitors are dynamically active 
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A. APPENDIX: Physical Basis of the Lorenz Model 



Mass conservation. Conservation of mass is enforced by use of a stream fu nction (ILandau &: Lifshitz 



( 19591 ). §9); the simplest solution is a two dimensional, cylindrical "roll" ( Chandrasekhai (1961) 



p. 44). The flow is assumed to be subsonic. 

Momentum conservation. The Navier-Stokes equation may be written as 

Du/Dt = i/V^u - gAp/p, (Al) 

in the low Mach number limit; the last term is the buoyant acceleration. The flow executes a circle 
of radius £/2; see Figure [3l The density fluctuation is related to the temperature fluctuation by 
Ap/po = /3t(T - To)/To, where /3t = -{d\np/d\iiT)p. 

We separate the variable u by considering the flow to be a mass flux which is constant around 
the ring at any given moment, represented by an average speed u which is a function of time onlj0. 
The hydrostatic background in the convective region has an entropy S which is constant. Then, 

IdP dW , 

--r = -i- = -ff' (A2) 

p dr dr 

where W = E + PV is the enthalpy per unit mass, and W = CpT, and Cp = -^^^IZY . Now 
J" gdr = gAr if g is constant. Lorenz assumed a linear temperature decrease with height, which 
corresponds to constant gravitational acceleration g. Using a height z = {£/2) cos (j), 

Wicp) = W{0) + g(i/2) cos (p, (A3) 

or 

T{cP)=To + ^^coscp. (A4) 

We denote Ti = g{l/2)/Cp, so 

Te = Tq + Ti cos (f), (A5) 

which corresponds to the environmental temperature. For Ti > g/Cp the system is convectively 
unstable, while for Ti < Tq, the background temperature can never be negative. We represent the 
temperature by 

T = To + T2 cos + T3 sin 0. (A6) 

Viscous damping. The viscous damping term may be approximated by i/V^u — >• —Fu = 
—u{2/ l)'^u. Here the constant F is the inverse of the viscous dissipation time scale T^is- 
The buoyant acceleration in the vertical direction is 

B, = -^i6p/p) ^5/3T(2T3/ro)sin2<^. (A7) 

Only a temperature (buoyancy) difference in the horizontal direction {(f) = tt/2) gives a net torque 
to turn the convective roll. 



Variations in density due to change in hydrostatic stratification could be included as well as the implied changes in 
velocity and cross-sectional area; see iTrittonI (|l988h . p. 188-196). 
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With damping due to a linear viscosity term, 

du/dt - 



ru+^{T-TE)sm. 



Integrating Eq. IA8I over a complete cycle in < (f) < 2tt, givea^^l 



du/dt 



-Tu + 



2n 



(A8) 



(A9) 



The temperature terms in cos (p integrate to zero because of the sin cj) factor in the buoyancy term. 
Thus we have just two terms in Eq. IA91 a sink from the viscous damping and a source from 
buoyancy. 

In the turbulent case, we might identify this rate of dissipation as an integral over the turbulent 
region, so it is transformed into a global quantity, ek/p = u'^/i- This is the deceleration times 
the velocity, giving a deceleration of —u\u\/£, so that i/V^u — )• —u\u\/£. The absolute value \u\ 
is used here because the characteristic time is t^s = ^/l^^l) and the deceleration is thus —u/rms- 
The length scale is the depth of the convective shell (i.e., £). T he convective speed u is constant in 
space, so we may use a nonlinear damping term, as implied bv iKolmogorovl (|l94ll ). 



2Tn 



(AlO) 



du/dt = —u\u\/£ + 

Energy conservation. At constant pressure, the first law of thermodynamics simplifies to 

dW/dt = CpdT/dt = ■¥ + (All) 



where e is the net heating-cooling rate from nuclear burning and neutrino emission, Cp the specific 
heat at constant pressure, T the temperature, p the mass (nucleon) density, F the flux of radiative 
energy, and ek the volumetric heating by turbulence. Ignoring e and £k, 



dT/dt + u • VT 



-{-V-F)/Cp. 
P 



The divergence of fiux of radiative luminosity is 



1 



-V • F 



1 



V-[-(c/3pK)VaT^], 



so that 



dT/dt + u • VT = z^tV^T, 



(A12) 



(A13) 



(A14) 



where vt = {^acT'^ /2>pk) / {pC p) is the thermal conductivity for radiative diffusion, for constant k 
and Cp. 

The background temperature Te is chosen such that V'^Te = 0. The radiative diffusion term is 



i/tV^T « K{Te-T), 



(A15) 



^The two terms constant in space (4>) generate a factor of 2it while the integral over sin'^ 4> gives a factor of tt. 
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where K{Te — T) is the extra radiative diffusion due to the deviation of temperature from Te, and 
K = VTi^/i)"^ is the inverse of the radiation diffusion time scale r,., 



ad- 



The classical Lorenz equations. Fohowing iLorenz we correct for an aspect ratio a by a 



factor 6 = 4/(1 + a^) which deals with the excess in vertical over horizontal heat conduction {bK 
versus K) . Substituting for T and Te in Equation IA141 we have 

dT2 .^dTs . 

—— cos m H — sm (p 

dt ^ dt ^ 

2uT2 . . , 2uT3 , 
— sm (p H — cos (p 

= bK{Ti-T2) cos (j) - KT3 sin (p. (A16) 

Coefficients of orthogonal functions must separately sum to zero, so 

dT3/dt-2uT2/e = -KT3 (A17) 
dT2/dt + 2uT3/i = bK{Ti-T2). (A18) 

We introduce a potential temperature 

n = Ti-T2, (A19) 

and eliminate T2. If Ti = gi/2Cp is independent of time, 

dT3/dt= -KT3 + 2uTi/i-2uTi/£ (A20) 
dT^/dt = -bKTi + 2uT3/e. (A21) 

We define dimensionless variables 

T = t K 

X = u {2 /IK) 

Y = T3 {g(3T/£rKTo) 

Z = {gl3T/irKTo), (A22) 

and use Equations I A91 lA20l and lA2ll to get the Lorenz equations in their usual form: Equations [T2t 
[la andm 

dX/dr = -aX + aY (A23) 

dY/dT = -XZ + rX-Y (A24) 

dZ/dr = XY - bZ, (A25) 

where r is a time in thermal diffusion units, a = T/K is the effective Prandtl number, and r = 
{g PtTi / iV KTq) = [gf^TTi/ lT'^TQ)a is the ratio of the Rayleigh number to its critical value for onset 
of convection. The Prandtl number is the ratio of coefficients of the viscous dissipation term to the 
thermal mixing term, a = v/vt = T/K = T^ad/'Tvisi where Tms = l/E is the viscous damping time. 
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B. APPENDIX: What Should The Prandtl Number Be? 



The original Lorenz parameters appear to be a fair choice to represent the flow in 3D simulations 
(see Figure [1]) . This is consistent with an effective Prandtl number for the numerical simulation of 
at ~ 10; however, the severe reduction in degrees of freedom from the simulations to the Lorenz 
model warns against taking the numerical value too literally. We note f or comp a rison that water 
has (T ~ 6 and air has a ~ 0.7. Apparently Lorenz felt he was lucky (jLorenj (|l993l l. p. 137); 



he took a suggested value (cr = 10) which gave chaotic behavior instead of a lower value, actually 
appropriate for air, for which his equations give stable rolls. As a mathematical example this 
quantitative difference in the Prandtl number is not significant, but for physical applications, it is. 

B.l. The Microscopic Value 

Using the ratio of microscopic thermal diffusion to vis cous time scale, iHansen &: Kawaler (Il994l ) 



(page 178 and 185) suggest th;t a « IG'^. S) finds the critical value for r for msta^y 



of steady convection to be Tc = a(a + 6 + 3) /(a — 6 — 1). For a < b + 1 = 2.666, steady convection 



is always stable, so that the [Hansen &: Kawaleii (|1994| ) value would never give turbulence. If 



a > 6 + 1, steady convection is unstable for sufficiently large Rayleigh numbers. This precise value 
for instability is a characteristic of the canonical Lorenz model, and is affected by the particular 
choice of dissipation function (see Appendix A). 

B.2. The Simulation Value 

In the numerical simulations, the effective Prandtl number is dominated by the turbulent cascade, 
which gives mixing of both momentum and heat at a rate determined by the largest eddy size. Thus, 
in numerical simulations, the turbulent flow deflnes its own effective value of this parameter a ^ at, 
and whatever value at has, it is clearly above the threshold for instability for the system we have 
numerically simulated. Numerical simulations and experiment suggest, for developed turbulence at 
high Reynolds numbers, at ~ 0.7, a value typical of many common gases. 

We note that in the "Reynolds analogy" ( Monin &: Yaglom ( 197ll ). p. 341), Osborne Reynolds 



argued that the mechanisms for transport of heat and momentum were essentially the same in a 
turbulent medium, so that the effective turbulent Prandtl number should be of order unity {at w 1). 

B.3. A Cascade Estimate 

Suppose we think of the Prandtl number as the relative strength of the process which makes the 
velocity field isotropic to that which converts the kinetic energy into heat. These occur at different 
ends of the turbulent cascade. To better understand what this might mean, consider a =(time 
to change direction) /(time to heat). We approximate the time to change direction by the time to 
halve the kinetic energy, {Ijlvrms)- At any length scale A the turbulent cascade has a transfer rate 
for kinetic energy oiv\l\ = v"^ jl. This implies a time spent at each length scale of t\ = X/v\. If 
each level in the cascade is smaller on average in length scale by a factor A(n + 1)/A(n) = /, the 
total time for the cascade is {i/v){l + f^/^ + • • •) = (i/v) /(I — /^^^), where the geometric series has 
been used for the summation (see also iFrisch ( 19951 ). §7.8). This gives an effective Prandtl number 



at ?s 2/(1 — f"^^^). Some representative values: for / = 1/e, this gives a ^ 4.1, for / = 1/2 this 
gives a 5, and / = 1/V2 gives a ~ 10. This argument may tend to overestimate the Prandtl 
number, so that at would be smaller than the historical value chosen by Lorenz (cr = 10). 
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B.4. Renormalization Group 



The separation in size of the large scale eddies with those at the dissipation scale suggests that 
this coupling might be approximated in some ingenious way. The microscopic viscosity might be 
considered a "bare" value, to be renor malized t o a "d ressed" value, in analogy to the field-theoretic 
treatment of interacting particles. iKadanofi ( 19661 ) proposed the idea of reducing the size of a 
system a step at a time by groupi ng neighborin g entities (in this case molecules) and treating 
each group as a single interaction. Wilson ( 197d ) has successfully implemented the general idea 
of coarse-graining, or "weeding out the sin all scales" . These general ideas may be applied to the 
turbulent cascade. Yakhot Qrszag ( 19861 ) use renormalization-group (RG) analysis of turbulence 
with some success, and in particul a r, est imate the effective Prandtl number for turbulent flow to 
be at = 0.7179; see also Kraichnan ( 198?! ). 

However there is some debate: in their review, Smith &: Woodrufl ( 19981 ) warn "...the RG 
method . . . leads to suggestive results when applied to turbulence. . . However, its application to 
turbulence cannot yet be called a major success, owing to th e uncontrolled a pproximations currently 
required to implement it." A similar sentiment is found in iMcCombI ( 20041 ) . p. 290. 



B.5. A Perspective 

We have captured the behavior of the pulses in the simulations, by a Lorenz model using the 
parameters that Lorenz used. In this sense, these values are relevant to our problem, although it is 
unclear from first principles what the precise value of the effective turbulent Prandtl number should 
be. The form for dissipation that Lorenz used is not the sam e as implied by Kolmogorov, so that the 
threshold for instability changes (|Arnett fc Meakinllioilc] ). We may interpret the Prandtl number 
a and the Rayleigh number r in terms of a renormalization in which the existence of turbulence 
implies effective Prandtl and Rayleigh numbers for the convective cell. The actual values of the 
microscopic viscosity and thermal conductivity have little feedback on the behavior of the largest 
eddies; see Section [2.11 Magnetic fields in real stars may affect the effective Prandtl number further, 
which in the classical Lorenz model affects in turn the development of instability in the convective 
roll. 

For problems which are insensitive to the details of the small scale flow, the values of the 
microscopic Prandtl do not matter; the turbulent system bootstraps to an effective dissipation 
for which the rules of Kolmogorov hold. This approximation applies to hydrostatic (and mildly 
dynamic) stellar evolution, in which the burning times are long compared to sound travel times. 
For these problems, ILES is appropriate. There are notable exceptions, such as a flame front 
in a medium of unmixed fuel and ash (SNIa progenitor models), for which the small scales are 
important, and direct numerical simulations (DNS) of the small scales is necessary. 

Finally, we recall how the Lorenz model was constructed: a low order mode was chosen, a 
Rayleigh number was chosen just above the onset of instability for that mode, and a Prandtl number 
was chosen which gave int eresting behavior. While different parameter choices are mathematically 
interesting ( Sparrow 19821 ). their physical relevance must be re-evaluated. The Lorenz equations 
are a spartan subset of the fluid equations which contain the germ of chaos; it is probably better 
to use them clS 3j guide rather than a gospel. 
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